Large  Deformation  Dynamic  Three-Dimensional  Coupled 
Finite  Element  Analysis  of  Soft  Biological  Tissues  Treated 

as  Biphasic  Porous  Media 


by  RA  Regueiro,  B  Zhang,  and  SL  Wozniak 


ARL-RP-0512 


November  2014 


Reprinted  from  CMES:  Computer  Modeling  in  Engineering  and  Sciences.  2014;98(l):l-39. 


Approved  for  public  release;  distribution  is  unlimited. 


NOTICES 

Disclaimers 

The  findings  in  this  report  are  not  to  be  construed  as  an  official  Department  of  the  Army  position  unless 
so  designated  by  other  authorized  documents. 

Citation  of  manufacturer’s  or  trade  names  does  not  constitute  an  official  endorsement  or  approval  of  the 
use  thereof 


Destroy  this  report  when  it  is  no  longer  needed.  Do  not  return  it  to  the  originator. 


Army  Research  Laboratory 

Aberdeen  Proving  Ground,  MD  21005-5069 


ARL-RP-0512 


November  2014 


Large  Deformation  Dynamic  Three-Dimensional  Coupled 
Finite  Element  Analysis  of  Soft  Biological  Tissues  Treated 

as  Biphasic  Porous  Media 


RA  Regueiro  and  B  Zhang 
University  of  Colorado 

SL  Wozniak 

Bowhead  Science  and  Technology,  LLC 


Reprinted  from  CMES:  Computer  Modeling  in  Engineering  and  Sciences.  2014;98(l):l-39. 


Approved  for  public  release;  distribution  is  unlimited. 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering 
and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  the  burden,  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson 
Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to 
comply  with  a  collection  of  information  if  it  does  not  display  a  currently  valid  OMB  control  number. 

PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  is  unlimited. 

13.  SUPPLEMENTARY  NOTES 

Reprinted  from  CMES:  Computer  Modeling  in  Engineering  and  Sciences.  2014;98(1):  1—39. 

14.  ABSTRACT 

The  report  presents  three-dimensional,  large  deformation,  coupled  finite  element  analysis  (FEA)  of  dynamic  loading  on  soft 
biological  tissues  treated  as  biphasic  (solid-fluid)  porous  media.  An  overview  is  presented  of  the  biphasic  solid-fluid  mixture 
theory  at  finite  strain,  including  inertia  terms.  The  solid  skeleton  is  modeled  as  an  isotropic,  compressible,  hyperelastic  material. 
FEA  simulations  include:  (1)  compressive  uniaxial  strain  loading  on  a  column  of  lung  parenchyma  with  either  pore  air  or  water 
fluid,  (2)  out-of-plane  pressure  loading  on  a  thin  slab  of  lung  parenchyma  with  either  pore  air  or  water  fluid,  and  (3)  pressure 
loading  on  a  l/8th  symmetry  vertebral  disc  (nucleus  and  annulus)  with  pore  water.  For  the  simulations,  mixed  formulation 
Q27P8  and  stabilized  Q8P8  finite  elements  are  compared  (“Q”  indicates  the  number  of  solid  skeleton  displacement  nodes,  and 
“P”  the  number  of  pore  fluid  pressure  nodes).  The  FEA  results  demonstrate  the  interplay  of  dynamics  (wave  propagation 
through  solid  skeleton  and  pore  fluid),  large  deformations,  effective  stress  and  pore  fluid  pressure  coupling,  compressibility  and 
viscosity  of  pore  fluid,  and  three-dimensional  effects  for  soft  biological  tissues  treated  as  biphasic  porous  media. 


15.  SUBJECT  TERMS 

soft  biological  tissues,  biphasic  mixture  theory,  dynamics,  large  deformations,  coupled  three-dimensional  finite  element  analysis 

17.  LIMITATION  18.  NUMBER 
OF  ABSTRACT  OF  PAGES 

UU  44 


19a.  NAME  OF  RESPONSIBLE  PERSON 

SL  Wozniak _ 

19b.  TELEPHONE  NUMBER  ( Include  area  code) 

410-306-2981 


16.  SECURITY  CLASSIFICATION  OF: 


a.  REPORT 

b.  ABSTRACT 

c.  THIS  PAGE 

Unclassified 

Unclassified 

Unclassified 

1 .  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE 

November  2014  Reprint 

4.  TITLE  AND  SUBTITLE 

Large  Deformation  Dynamic  Three-Dimensional  Coupled  Finite  Element  Analysis 
of  Soft  Biological  Tissues  Treated  as  Biphasic  Porous  Media 


6.  AUTHOR(S) 

RA  Regueiro,  B  Zhang,  and  SL  Wozniak 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

University  of  Colorado  Boulder 
1 1 1 1  Engineering  Dr 
Boulder,  CO  80309 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

US  Army  Research  Office 
POBox  12211 

Research  Triangle  Park,  NC  27709 


3.  DATES  COVERED  (From  -  To) 

7  June  2011-20  September  2014 

5a.  CONTRACT  NUMBER 

W91  INF-01  l-D-0001 
W91 1QX-09-C-0057 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

ARL-RP-0512 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

ARO 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

TCN  11-024 


11 


Standard  Form  298  (Rev.  8/98) 
Prescribed  by  ANSI  Std.  Z39.18 


Copyright  ©  2014  Tech  Science  Press 


CMES,  vol.98,  no.l,  pp.1-39,  2014 


Large  Deformation  Dynamic  Three-Dimensional  Coupled 
Finite  Element  Analysis  of  Soft  Biological  Tissues  Treated 
as  Biphasic  Porous  Media 


R.A.  Regueiro1,2,  B.  Zhang1 2,  S.L.  Wozniak3 


Abstract:  The  paper  presents  three-dimensional,  large  deformation,  coupled  fi¬ 

nite  element  analysis  (FEA)  of  dynamic  loading  on  soft  biological  tissues  treated  as 
biphasic  (solid-fluid)  porous  media.  An  overview  is  presented  of  the  biphasic  solid- 
fluid  mixture  theory  at  finite  strain,  including  inertia  terms.  The  solid  skeleton  is 
modeled  as  an  isotropic,  compressible,  hyperelastic  material.  FEA  simulations  in¬ 
clude:  (1)  compressive  uniaxial  strain  loading  on  a  column  of  lung  parenchyma 
with  either  pore  air  or  water  fluid,  (2)  out-of-plane  pressure  loading  on  a  thin  slab 
of  lung  parenchyma  with  either  pore  air  or  water  fluid,  and  (3)  pressure  loading 
on  a  l/8th  symmetry  vertebral  disc  (nucleus  and  annulus)  with  pore  water.  For 
the  simulations,  mixed  formulation  Q27P8  and  stabilized  Q8P8  finite  elements  are 
compared  (“Q”  indicates  the  number  of  solid  skeleton  displacement  nodes,  and 
“P”  the  number  of  pore  fluid  pressure  nodes).  The  FEA  results  demonstrate  the 
interplay  of  dynamics  (wave  propagation  through  solid  skeleton  and  pore  fluid), 
large  deformations,  effective  stress  and  pore  fluid  pressure  coupling,  compressibil¬ 
ity  and  viscosity  of  pore  fluid,  and  three-dimensional  effects  for  soft  biological 
tissues  treated  as  biphasic  porous  media. 

Keywords:  soft  biological  tissues;  biphasic  mixture  theory;  dynamics;  large  de¬ 
formations;  coupled  three-dimensional  finite  element  analysis 

1  Introduction 

It  is  well-known  that  soft  biological  tissues  are  multiphase  (oftentimes  treated  as 
a  biphasic  mixture  of  solid  and  fluid  phases  (Holmes,  1986;  Suh,  Spilker,  and 
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Holmes,  1991;  Almeida  and  Spilker,  1998;  Levenston,  Frank,  and  Grodzinsky, 
1998))  and  can  undergo  large  deformations.  Few  researchers  have  considered  in¬ 
ertia  terms  in  the  biphasic  theory  in  studying  wave  propagation  through  soft  bio¬ 
logical  tissues  (Zhu  and  Suh,  2000,  2001).  The  flow  of  pore  fluid  relative  to  solid 
skeleton  deformation  (such  as  squeezing  a  saturated  sponge)  must  be  handled  prop¬ 
erly  within  a  finite  strain  context,  as  well  as  considering  inertia  terms  for  dynamic 
loading  when  necessary,  and  implemented  properly  within  a  mixed  Lagrangian  fi¬ 
nite  element  (FE)  formulation,  or  other  large  deformation  numerical  method.  For 
typical  physiological  dynamic  loading  such  as  encountered  during  normal  athlet¬ 
ic  activities  (e.g.,  running,  jumping),  the  relative  acceleration  of  the  fluid  phase 
af  with  respect  to  solid  phase  as  may  be  approximated  as  zero:  af  =  af  —  as  ^  0. 
At  higher  strain  rates,  however,  such  as  encountered  during  shock  loading  and  head 
impact  (e.g.,  leading  to  traumatic  brain  injury  (TBI)),  the  relative  acceleration  of  the 
fluid  phase  with  respect  to  solid  phase  may  not  be  zero  (af  ^  0),  requiring  reformu¬ 
lation  of  the  balance  equations  originally  formulated  for  lower  rate  loading  (such 
as  normal  athletic  activities).  Such  extension  is  discussed  in  the  paper,  but  all  FE 
results  currently  assume  af  «  0.  The  mixed  FE  formulation  and  three-dimensional 
(3D)  Q27P8  hexahedral  FE  implementation  (27  solid  skeleton  displacement  nodes, 
8  pore  fluid  pressure  nodes,  Fig.4)  leads  to  a  stable  finite  element  method  even  for 
undrained  loading  conditions,  such  as  those  encountered  at  the  initial  transient  of 
dynamic  loading.  A  stabilized  mixed  formulation  (Brezzi  and  Pitkaranta,  1984; 
Truty  and  Zimmermann,  2006;  White  and  Borja,  2008;  Sun,  Ostien,  and  Salinger, 
2013)  Q8P8  hexahedral  element  is  also  implemented  within  the  coupled  dynamics 
framework,  and  results  are  compared  with  the  Q27P8  element.  The  formulation  for 
biphasic  mixture  theory  at  finite  strain  naturally  calculates  the  build  up  of  pore  flu¬ 
id  pressure,  and  thus  properly  calculates,  through  the  effective  stress  principle,  the 
change  in  solid  skeleton  stress  (i.e.,  the  “effective”  stress)  over  time,  when  a  bipha¬ 
sic  soft  tissue  is  subjected  to  dynamic  loading.  Also,  after  the  initial  transient,  the 
variation  of  solid  skeleton  stresses  will  be  naturally  calculated  as  the  fluid  phase 
pressure  dissipates  over  time.  This  is  important  for  developing  physiologically- 
relevant  degradation/damage  hyperelastic,  anisotropic  constitutive  models  for  soft 
biological  tissues  (Pena,  2011;  Balzani,  Brinkhues,  and  Holzapfel,  2012)  within 
the  context  of  biphasic  mixture  theory.  The  3D  FE  implementation  is  conducted 
in  Tahoe  (tahoe .  sourcef  orge .  net),  an  opensource  C++  FE  code,  with  more  de¬ 
tails  provided  in  Ebrahimi  (2007);  Regueiro  and  Ebrahimi  (2010),  and  formulation 
details  in  Li,  Borja,  and  Regueiro  (2004). 

An  outline  of  the  remainder  of  the  paper  is  as  follows:  Section  2  presents  a  brief 
overview  of  the  theory  of  biphasic  solid-fluid  mixtures  at  finite  strain,  including 
kinematics,  balance  equations  (linear  momentum,  and  mass),  and  thermodynamics 
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for  constitutive  equation  forms.  Section  3  also  presents  the  stabilization  term  on  the 
balance  of  mass  for  stabilized  Q8P8  implementation.  Section  4  presents  numerical 
examples  to  test  the  performance  of  the  Q27P8  versus  Q8P8  elements,  and  also  to 
study  the  importance  of  including  inertia  terms  for  simulating  dynamic  loading  of 
soft  biological  tissues  treated  as  porous  media.  Section  5  summarizes  the  results, 
conclusions,  and  future  work. 

Index  notation  will  be  used  wherever  needed  to  clarify  the  presentation.  Carte¬ 
sian  coordinates  are  assumed,  so  all  indices  are  subscripts,  and  the  partial  spatial 
derivative  is  the  same  as  covariant  spatial  derivative  (Eringen,  1962).  Some  sym¬ 
bolic/direct  notation  is  also  given,  such  that  ( FT  •  F)jj  =  F^Fu,  where  F  is  the 
deformation  gradient.  Boldface  denotes  a  tensor  or  vector.  Subscript  (#)5;  implies 
a  partial  spatial  derivative.  Lowercase  subscript  i  denotes  a  leg  of  the  tensor  in  the 
current  configuration  38,  and  capital  subscript  7(s)  denotes  a  leg  of  the  tensor  in  the 

def 

reference  configuration  38\  of  the  solid  skeleton.  Superposed  dot  (□)  =  Ds(U)/Dt 
denotes  material  time  derivative  with  respect  to  the  solid  skeleton  motion.  The 

dcf 

symbol  =  implies  a  definition. 

2  Biphasic  (solid-fluid)  mixture  theory  at  finite  strain:  overview  of  theory 
and  3D  FE  implementation 

The  background  for  the  theory  of  porous  media  at  finite  strain  may  be  found  pri¬ 
marily  in  Bowen  (1980,  1982);  Coussy  (2004);  de  Boer  (2005),  and  originally  in 
Truesdell  and  Toupin  (1960).  For  other  details  on  nonlinear  solid  mechanics,  refer 
to  Holzapfel  (2000)  and  references  therein.  We  follow  the  notation  of  Holzapfel 
(2000);  de  Boer  (2005),  and  to  some  extent  also  Bowen  (1980,  1982). 

2.1  Concept  of  volume  fraction  and  mixture  theory 

The  concept  of  volume  fraction  is  illustrated  in  Fig.  1  for  the  lung  parenchyma  (alve¬ 
olar  tissue  for  solid  skeleton).  The  volume  fractions  na  for  a  biphasic  mixture  (solid 
(s)  and  fluid  (f))  relate  “real”  quantities  with  respect  to  the  differential  volume  dva 
of  constituent  a  in  the  current  configuration,  versus  the  smeared  quantity  over  the 
total  differential  volume  dv,  where  na  —  dva/dv,  or  dva  =  nadv  .  For  example,  the 
partial  mass  density  of  the  a  constituent  is  calculated  as  pa  =  paRna  (see  Fig.l), 
where  paR  is  the  real  mass  density  of  constituent  a.  Similarly,  mass  of  constituent 
a,  ma,  over  the  total  body  38  can  be  defined  (see  Fig.l). 

2.2  Motion  and  kinematics ,  material  time  derivative 

The  kinematics  of  a  biphasic  solid-fluid  mixture  theory  are  shown  in  Fig. 2.  The 
vector  x  is  the  spatial  position  vector,  which  is  simultaneously  occupied  by  all 
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*theory  of  porous  media  assumes 
control  space  is  that  of  the  solid 
phase  38  =f  38s 


na(x,t)  =  dva/dv 

^ na  =  1  ,  nf  +  ns  =  1 
a 

dv  =  Y,dvoc 
a 

•  na  =  volume  fraction  of  constituent  a  in  dv  C 
38,  where  38  =  38s 

•  dva  =  differential  volume  of  constituent  a  in 
dv 

ma=  [  paRdva=  [  paRnadv  =  [  padv 

paR  =  dma/dva 
pa  =  dma/dv  = 

•  =  differential  mass  of  constituent  a 

•  paR  =  real  mass  density  of  constituent  a 

•  pa  =  partial  mass  density  of  constituent  a 


Figure  1:  Concept  of  volume  fraction  for  biphasic  (solid(s)-fluid(f))  mixture  theory, 
showing  solid  skeleton  composed  of  alveolar  tissue,  and  definitions  of  mass  and 
density. 


constituent  material  points  Xs,Xf  of  the  mixture  (homogenized,  or  smeared),  such 
that  x  =  =  £s(Xs,f),  where  the  material  point  of  the  solid  skeleton  Xs  is 

mapped  from  the  reference  position  Xs  to  the  current  position  x  through  mapping 
Xs  (and  similarly  for  the  material  point  of  pore  fluid  Xf  which  maps  through  Xf, 
Fig. 2).  We  define  the  inverse  map  Xa  =  ZaH***)  ( xi{a )  =  assuming 

smoothly  differentiable  fields.  The  deformation  gradient  and  its  inverse  are  written 
for  each  phase  a  as, 


Fa  = 


d Xa 

dXa  ’ 


dXq 

dx 


Fii(oc) 


^Xi(a) 
dXI{a)  ’ 


1  _  dXlja) 
II ^  dxi 


(1) 


Likewise,  the  volumetric  deformation  of  a  solid- fluid  mixture,  that  is  smeared  in  the 
current  configuration  at  spatial  position  vector  x,  is  shown  in  Fig. 3.  The  differential 
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Figure  2:  Kinematics  of  a  biphasic  (solid-fluid)  mixture  theory,  showing  solid 
skeleton  composed  of  alveolar  tissue.  The  continuum  assumption  of  mixture  theory 
is  evident  in  the  assumption  that  solid  (s)  and  fluid  (f)  constituents  coexist  at  the 
current  position  x ,  although  their  velocities  vs  and  Vf  may  be  different;  i.e.,  Vf  ^  vs, 
in  general. 


volumes  dV f  and  dVs  in  their  respective  reference  configurations  and  ^q,  both 
map  to  the  same  differential  volume  dv  in  the  current  configuration  £%,  through 
their  deformation  gradients  F f  and  Fs. 

The  Jacobian  of  deformation  for  the  two  constituents  is  written  as, 


Js  =  det.Fs  >  0  ;  J(  =  detFf  >  0 

(2) 

dv  =  JsdVs  =  JfdVf 

(3) 

dva  =  nadv  =  naJadVa 

(4) 

dVf  C  ,  dV s  C 

(5) 

where  we  will  typically  drop  the  s  superscripts  and  subscripts  because  the  theory 
of  porous  media  assumes  we  follow  the  motion  of  the  solid  skeleton. 
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Figure  3:  Volumetric  deformation  of  solid  and  fluid  constituents  in  a  biphasic  mix¬ 
ture  (solid  skeleton  composed  of  alveolar  tissue  of  the  lung  parenchyma). 


We  convert  all  material  time  derivatives  with  respect  to  the  solid  phase  motion 

•  ciof 

((□)  =  DS(D) /Dt),  such  that  for  phase  a,  the  material  time  derivative  is, 


£>«(□)  _  DS(D)  <?(□) 

i  ^  ‘  Va 


Dt 

Va  =  Va 

DS(D) 


Dt 

-Vs 


dx 


(6) 

(7) 

(8) 


Dt  dt  dx 
where  va  is  the  relative  velocity  vector  of  the  a  phase  with  respect  to  the  solid  (s) 
phase  motion.  The  material  time  derivative  will  be  used  in  deriving  the  balance 
equations  in  the  following  sections. 


2.3  Balance  of  mass  (spatial  and  material  descriptions) 

For  the  balance  of  mass  of  the  mixture,  we  write  separately  the  balance  of  mass  of 
each  constituent,  solid  and  fluid,  expressing  all  material  time  derivatives  in  terms 
of  the  solid  skeleton  motion,  and  then  add  the  two  equations  together  to  obtain  the 
balance  of  mass  of  the  mixture.  The  total  mass  of  constituent  a  in  is  written  as 
(cf.  Fig.l) 

ma=  [  padv  =  f  paJadVa 


(9) 
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Taking  the  material  time  derivative  of  this  spatial  field  ma  with  respect  to  the  motion 
of  constituent  a,  we  can  express  the  balance  of  mass  of  constituent  a  as 


Dama 

Dt 


=  /  (^+p<Mivv„U=  /  ^ 

JSS\  Dt  )  J ® 


(10) 


where  y®  is  the  mass  supply  rate  for  constituent  a.  If  we  assume  compressibility 
of  constituent  a,  assuming  temperature  is  constant,  we  may  write  (Li,  Borja,  and 
Regueiro,  2004) 


1  DapaR  _  1  Dapa 
pa*  Dt  ~  Ka  Dt 


(ID 


where  pa  is  the  mean  Cauchy  stress  of  constituent  a,  and  Ka  is  the  bulk  modulus. 
Then  balance  of  mass  for  constituent  a  becomes: 


Dan° 


Dt 


rPWjpa 

Ka  Dt 


+  nadivv«  = 


yx 

naR 


(12) 


We  assume  that  the  solid  constituent  is  nearly  incompressible,  such  that  ns/Ks  — > 
0,  and  also  that  there  is  no  supply  of  solid  mass,  such  that  7s  =  0.  Using  these 
assumptions,  we  can  solve  for  the  volume  fraction  of  the  solid  phase  as  ns  =  ns0/Js 
in  closed  form  from  the  balance  of  mass  equation  for  the  solid  phase  by  itself  (Li, 
Borja,  and  Regueiro,  2004).  The  volume  fraction  of  fluid  is  then  solved  as  /;'  = 
1  —  ns.  The  simplified  form  of  the  balance  of  mass  of  the  biphasic  (solid-fluid) 
mixture  (to  solve  for  Cauchy  pore  fluid  pressure  pi)  then  results  as,  when  adding 
Eq.(12)  for  a  =  f,  s, 


n{  Dsp{  1  dpi  .  f_  .  ,  f_  .  / 

K,  or + d,vv- + K,  li '  (" Vf) + d,v(" v,)  =  7* 


(13) 


We  can  map  Eq.(13)  back  to  the  reference  configuration  of  the  solid  phase 
obtain  the  following, 


to 


V  D>r  ,  DSJS  (  Js  dpt 


+ 


K(  Dt  Dt  KfdX 


F 


-1 


(nfVf)  +JS 


dX. , 


.  p-T  Q4) 

•  r  s  pfR 


For  a  Total  Lagrangian  FE  implementation,  this  is  the  equation  from  which  we 
derive  our  variational  equation  for  the  weak  form. 
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2.4  Balance  of  linear  momentum  ( material  and  spatial  descriptions ) 

Here,  the  balance  of  linear  momentum  of  the  biphasic  solid-fluid  mixture  is  pre¬ 
sented.  Carrying  out  the  material  time  derivative  of  the  linear  momentum  with 
respect  to  the  a  phase  motion,  applying  the  balance  of  mass  of  constituent  a,  and 
using  the  divergence  theorem  on  the  traction  term  ta,  we  can  localize  the  integral 
to  obtain  the  balance  of  linear  momentum  for  phase  a  in  the  current  configuration 
as, 

di  \Ga  +  paba  +  ha  =  paaa  +  f^Va  (15) 


where  Ga  is  the  partial  Cauchy  stress,  such  that  the  total  Cauchy  stress  G  =  (7s  + 
<7f,  pa  is  the  partial  mass  density,  such  that  total  mass  density  p  =  ps  +  pf,  ba  is 
the  body  force  per  unit  mass  on  constituent  a  (we  will  assume  the  same  body  force 
per  unit  mass  for  each  constituent,  such  as  acceleration  of  gravity,  ba  =  g ),  ha  is 
the  interaction  body  force  from  all  other  constituents  on  constituent  a,  aa  is  the 
acceleration  vector,  and  is  the  mass  supply  momentum  (usually  negligible). 
We  note  that  the  internal  body  forces  due  to  drag  between  constituents  sum  to 
zero  (equal  and  opposite),  and  thus  do  not  affect  the  mixture  as  a  whole,  such  that 
hs  +  hf  =  0.  The  balance  of  angular  momentum  for  non-polar  constituents  states 
that  the  respective  partial  stresses  (and,  in  turn,  the  total  Cauchy  stress,  and  effective 
stress)  are  symmetric:  Ga  =  (Ga)T . 

It  is  now  relevant  to  discuss  a  principle  that  allows  us  to  distinguish  stress  acting  on 
the  solid  skeleton,  and  the  pressure  acting  on  the  pore  fluid  (assuming  the  fluid  is 
nearly  inviscid,  such  as  water).  We  apply  the  effective  stress  principle,  which  can 
be  credited  to  Terzaghi  (1943)  (pg!2)  for  saturated  condition  of  soils,  that  states1 


G  =  Gf  -  pf 


j^skel 


i 


(16) 


where  the  real  pore  fluid  pressure  pf  =  ^tr(af)  is  positive  in  compression,  and  the 
mean  effective  stress  is  positive  in  tension  p'  —  -Itr!  <7  ).  where  <7  is  the  “effective” 
Cauchy  stress,  or  the  stress  acting  on  the  solid  skeleton,  for  which  we  will  apply 

(r^skd  \ 

1  —  J  is  the  Biot  coefficient 


(Coussy,  2004),  and  Kskel  is  the  solid  skeleton  bulk  modulus.  For  soils  and  soft 
biological  tissues,  Kskel/Ks  0,  whereas  for  rocks  and  bone  Kskel/Ks  is  finite.  We 


1  The  application  of  the  effective  stress  principle  to  soft  biological  tissues  has  not  been  completely 
tested  to  date,  but  it  is  applied  here  for  theoretical  and  numerical  convenience.  This  is  a  topic  of 
further  research. 
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assume  a  nearly  inviscid  (no  shear  stress)  isotropic  fluid  (e.g.,  water),  where  then 


G 


f 


(17) 


and 


(18) 


Thus,  we  note  that  the  partial  solid  stress  Gs  is  not  equal  to  the  effective  stress  g' , 
Gs  ^  G ',  unless  the  pore  fluid  pressure  pf  =  0.  The  effective  stress  principle  is 
useful  for  introducing  constitutive  equations  for  the  solid  skeleton  separate  from 
the  pore  fluid. 

Starting  with  Eq.(15),  we  can  map  the  balance  of  linear  momentum  of  phase  a 
back  to  the  reference  configuration  of  the  solid  phase  The  partial  first  Piola- 
Kirchhoff  stress  a  with  respect  to  solid  phase  reference  configuration  is  written 
as, 

pf  =  jsoa  •  f;t  ,  J*(s)  =  (19) 

where  subscript  s  denotes  the  reference  configuration  38^  to  which  the  j  leg  of  G® 
is  mapped.  It  is  then  possible  to  arrive  at  the  balance  of  linear  momentum  in  the 
reference  configuration  of  the  solid  phase  (s)  for  phase  a  as, 

DIVs  Pf  +  p^ba  +  Jsha  =  p0«(s)aa  +  ^(s)v«  (20) 

Starting  with  Eq.(20),  we  can  write  each  balance  of  linear  momentum  equation  for 
solid  (s)  and  fluid  (f)  phases,  and  use  the  following  information  to  derive  the  balance 
of  linear  momentum  of  the  solid-fluid  mixture  (s)+(f)  in  the  reference  configuration 
of  the  solid  phase  &o- 


1.  total  Cauchy  stress,  and  first  Piola-Kirchhoff  stress  with  respect  to  &o- 

a  =  <7s  +  <Tf,  P,  =  P\  +  P{  (21) 


2.  effective  stress  equation  for  Cauchy  stress: 

/  t^skel N 

a  =  a’  -  pf  f  1  - 


K, , 


1 


(22) 
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3.  assume  solid  and  fluid  phase  accelerations  are  nearly  the  same  (for  now): 
df  as  =  a ,  which  may  be  appropriate  for  longer  period  motions  like  earth¬ 
quakes  and  athletic  activities,  but  likely  not  appropriate  for  high  impact  events 
experienced  during  car  crash,  or  blast  loading 

4.  assume  all  mass  supplies  are  negligible:  y01  —  0 

5.  assume  body  forces  per  unit  mass  are  only  due  to  gravity:  ba  =  g ,  where  g 
is  the  acceleration  vector  of  gravity 

The  resulting  balance  of  linear  momentum  of  the  solid-fluid  mixture  (s)+(f)  in  the 


reference  configuration  of  the  solid  phase  is  then  written  as, 

DFVs^s  +  Po(s)£  =  Po(s)a  (23) 

Ps  =  P's-JsP{BF;t  (24) 

Eventually,  for  the  Total  Lagrangian  finite  element  formulation,  we  will  drop  the 
s  designation  because  the  reference  configuration  will  always  be  that  of  the  solid 
skeleton  (phase),  such  that, 

DIVP  +  pog  =  poa  (25) 

P  =  P'  -JpfBF~T  (26) 


Remark  1:  balance  of  linear  momentum  (spatial  description)  for  as  ^  a f.  Let 

us  revisit  the  balance  of  linear  momentum  for  the  solid-fluid  mixture  in  the  current 
configuration  assuming  as  ^  df,  such  that, 

diver  +  pg  =  psas  +  pfaf  (27) 

The  question  becomes  how  to  handle  df.  One  way  is  to  write  the  balance  of  linear 
momentum  specifically  for  the  fluid  phase  as, 

di  \Gf  +  pfbf  +  hf  =  pfdf  (28) 

where  for  an  inviscid  fluid  phase,  <7f  =  —  nfpf  1,  we  can  express  the  interaction  fluid 
body  force  as, 

f  dnf  f  nf\  f_ 


(29) 
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where  the  superficial  fluid  velocity  vector  (nfv f)  (a.k.a.,  the  Darcy  velocity)  is  de¬ 
fined  constitutively  in  Eq.(41).  Then,  upon  substituting  Eq.(29)  into  Eq.(28),  we 
can  write  the  balance  of  linear  momentum  for  the  fluid  phase  as, 

pfaf  +  nf^  +  (j)  (”fi>f)  ~PfS  =  0  (30) 

Then,  Eq.(30),  upon  mapping  back  to  the  solid  skeleton  reference  configuration, 
including  Eqs.(14),(23),  we  have  three  coupled  balance  equations  to  solve  for  three 
unknown  fields:  solid  skeleton  displacement  us ,  fluid  phase  displacement  Hf,  and 
Cauchy  pore  fluid  pressure  pf.  A  similar  procedure  was  followed  for  small  strain 
by  Jeremic,  Cheng,  Taiebat,  and  Dafalias  (2008).  This  is  left  for  future  work. 

2.5  Thermodynamics  (first  and  second  laws ,  constitutive  equation  forms) 

Before  we  introduce  constitutive  equations,  we  briefly  present  the  thermodynamics 
for  a  biphasic  mixture.  Using  the  first  and  second  laws  of  thermodynamics,  and 
assuming  the  existence  of  a  Helmholtz  free  energy  per  unit  mass  of  the  a  phase 
i//a,  we  derive  the  Clausius-Duhem  inequality,  which  will  be  useful  for  defining 
constitutive  model  forms. 

Applying  the  material  time  derivative  to  the  total  internal  energy  of  phase  a,  using 
the  balance  of  mass  for  a,  divergence  theorem  on  the  traction  term,  and  balance  of 
linear  momentum  on  a,  and  localizing  the  integral,  we  can  derive  the  balance  of 
energy  of  phase  a  in  the  current  configuration  SB  as, 


r>a  a 

pa -£a:Ga  +  di  vqa  -  para 

\  .  stresspower  internal  heat  flux  internal  heat  supply  rate 

internal  energy  density  rate 

=  f  Qva-Va-e0^  -  (31) 

s - v« - y  interphase  power  phase  a  power 

mass  supply  power 


where  ea  is  the  internal  energy  per  unit  mass  of  a,  qa  is  the  heat  flux  vector,  ra 
is  the  heat  input  rate  per  unit  mass,  and  ea  is  the  power  density  supply  to  phase  a 
by  other  phases.  The  second  law  of  thermodynamics  for  phase  a  is  written  in  the 
current  configuration  SB  as, 


Dana 

9a,yaria  +  pa - —6a  —p 


a  ra 


i  dea 
~e^~dx 


■  qa  +  d\\qa  >  0 


Dt 


(32) 
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where  V)a  is  the  entropy  per  unit  mass,  and  0a  is  the  temperature  of  phase  a.  To 
obtain  the  Clausius-Duhem  inequality  in  the  current  configuration  for  phase  a,  we 
consider  the  existence  of  the  Helmholtz  free  energy  per  unit  mass  for  phase  a  in 
the  current  configuration  y/a  as, 

\l/a  =  ea  -Oaria  (33) 


Taking  the  material  time  derivative  of  Eq.(33),  using  the  first  law  in  Eq.(31),  and 
substituting  into  the  second  law  in  Eq.(32),  we  arrive  at  the  Clausius-Duhem  in¬ 
equality  for  phase  a  as, 


mass  supply  power 

Da  \i/a 

p  - — 

H  Dt 


~  ( Pari 


stresspower  interphase  power  phase  fx  power 

Daea  i  dea 


Dt 


9a  dx 


■q  >0 


(34) 


Helmholtz  free  energy  density  rate  temperature  rate  power  heat  flux  power 


We  assume  the  solid  constituent  is  nearly  incompressible  ( psR  is  constant)  and  the 
fluid  phase  is  compressible,  such  that  the  Helmholtz  free  energy  per  unit  mass  for 
the  solid  skeleton  (s)  in  and  fluid  phase  (f)  in  SS  are  written  as, 

p^s(Cs,0s)  ;  pW*,ef)  (35) 

where  Cs  =  F\  •  Fs  is  the  right  Cauchy-Green  tensor  for  the  solid  skeleton  de¬ 
formation.  We  assume  the  fluid  partial  stress  Gf  consists  only  of  a  pressure  term 
(inviscid,  such  as  water),  such  that  Gf  —  —pf  1,  where  the  partial  fluid  pressure  pf  is 
related  to  the  Cauchy  fluid  pressure  pf  through  pf  =  nfpf.  Recall  the  total  Cauchy 
stress  written  in  terms  of  the  partial  stresses  and  solid  skeleton  effective  stress  as, 

a  =  os  +  of  =  o'  -  pfBl  (36) 

We  write  the  Clausius-Duhem  inequality  for  each  phase  ( a  =s,f)  in  Eq.(34),  add 
them  together,  account  for  the  functional  forms  of  the  Helmholtz  free  energy  func¬ 
tions  in  their  respective  configurations  in  Eq.(35)  (^q  and  respectively),  and 
then  derive  the  Clausius-Duhem  inequality  for  the  solid-fluid  mixture  as  (all  terms 
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pushed  forward  to  the  current  configuration  and  localizing  the  integral), 


1  _/  X 

—S' - 


1  d[p^y\  DSCS 


Js  dCs 


+  [PfB^m 


d(pV)  1  DfP 


+  •  («f%) 


1  d(ps0Ys)  1  ,  s  s, 


V  W  „  1  KJ  f 

-^q-¥^xq^ 


<3(pV) 


+  (p  Y) 


0s  dx 


Applying  the  Coleman  and  Noll  (1963)  argument  for  independent  rate  processes 
and  then  for  Eq.(37)  to  be  satisfied  (and  hence,  the  2nd 


Dt  ’  Dt  ’  Dt  ’  Dt  ’  AWA  ^  »■ 

law),  the  following  constitutive  equations  must  hold: 


PfR  d(PY 


c,  _^d{PoVs)  _ptRd(pV) 

Sa~2~dcTtPi-Bif~d^~  (: 

oSt?s_  d(PlV)  f  f_  *(pV)  r 

P0T?  -  ^0^  ’ p  p  -  ^0^  (- 

The  remaining  terms  in  Eq.(37)  comprise  the  reduced  dissipation  inequality  as, 


Bn(  dpm 

f„f  <YV) 


f  f 

,  pY  =  - 


dPf  ,  f R(  .  f 


0s 


0f 


q  >0 


Furthermore,  we  assume  thermodynamic  conjugacy  through  proportionality  pa¬ 
rameters  (permeability  k ,  and  thermal  conductivities  ke s  and  ),  such  that  the 
following  constitutive  forms  hold  for  generalized  Darcy’s  law,  and  Fourier’s  law 
(assuming  local  thermal  equilibrium  such  that  mixture  temperature  0  =  0s  =  0f): 

generalized  Darcy’s  law: 

Y~  def  \*Pf  ,  Rr~  tf\l 


UCi  ,  /  t 

n  Vf  =  —  k{n 


(af  —  b) 


T]f5(^0)  l-(^)z 


where  x  is  the  intrinsic  permeability,  r/f  is  the  fluid  viscosity,  and  <5(nf)  is  the 
Cozeny-Karman  relation  for  porosity  dependence  of  k  (as  a  function  of  solid  skele¬ 
ton  volume  change  Js  as  nf  =  1  —  ns ,  ns  =  ns0/Js),  with  /i[}  the  initial  porosity 
(Coussy,  2004). 
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Fourier's  law : 


q  =  cf  +  qf  d=lf  d-2  ,  kemx  =  nskdS  +  nfkei 

OX 


(42) 


We  now  consider  specific  equations  for  the  Helmholtz  free  energy  functions: 

Fluid :  we  assume  the  following  form  for  the  Helmholtz  free  energy  function  for 
the  fluid: 


pY(p«0f)  =  lB(nfK{)(lnpfR)2+gf(ef) 


(43) 


where  gf  (0f)  is  a  temperature-dependent  term,  if  needed.  We  assume  homogeneous 
temperature,  isothermal  conditions  for  now.  Using  Eq.(38)2,  we  can  derive  the  con¬ 
stitutive  equation  for  the  real  mass  density  of  the  fluid  pfR  in  terms  of  the  Cauchy 
pore  fluid  pressure  pf  as, 


P 


fR 


Po^exp 


Pf  ~  Pm 
Kt 


(44) 


Solid  skeleton :  we  assume  the  following  form  for  the  Helmholtz  free  energy  func¬ 
tion  for  the  solid  skeleton  (neo-Hookean  compressible  isotropic  elasticity  (Ogden, 
1984)): 

pW(Cs,ds)  =  U(Js)  +  ^p(trCs-3)+gs0(Gs) ,  U(JS)  =  ^l(lnJs)2-p(lnJs)  (45) 

where  go(®S)  a  temperature-dependent  term.  Using  Eq.(38)i,  we  then  derive  the 
constitutive  equation  for  the  effective  Second  Piola-Kirchhoff  stress  as, 

5' =  [A0n7s)-p]Cs7'+pl  (46) 


By  assuming  functional  form  Pq1//s(Cs,  0s)  for  the  solid  skeleton,  we  can  later  gen¬ 
eralize  for  anisotropy  (fiber  directions)  of  soft  tissues  (Holzapfel  and  Gasser,  2000). 


3  Stabilized  Finite  Element  Implementation 

The  nonlinear  finite  element  formulation  and  implementation  (Newton-Raphson 
nonlinear  solution,  Newmark  time  integration)  is  discussed  in  Li,  Borja,  and  Regueiro 
(2004);  Ebrahimi  (2007);  Regueiro  and  Ebrahimi  (2010),  and  thus  details  are  not 
presented.  We  will  focus  on  the  stabilized  term  to  be  added  to  the  variational  equa¬ 
tion  for  the  balance  of  mass  of  the  solid-fluid  biphasic  mixture  as  discussed  in  the 
references. 
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3.1  Stabilized  term 

For  stabilizing  the  tri-linear  hexahedral  Q8P8  element  near  the  undrained  condition 
(upon  transient  loading,  with  low  permeability,  and  assuming  nearly  incompress¬ 
ible  fluid  phase  Kf  — »  oo)?  we  follow  the  approach  of  Truty  and  Zimmermann  (2006), 
which  is  based  on  the  method  of  Brezzi  and  Pitkaranta  (1984).  We  attempted  to  ap¬ 
ply  the  stabilization  procedure  of  White  and  Borja  (2008)  for  projecting  pore  fluid 
pressure  pf  by  integral-averaging  over  an  element  e ,  along  with  its  weighting  func¬ 
tion  value  r] ,  but  based  on  our  implementation,  it  was  ineffective  for  our  particular 
applications  of  soft  tissues  at  finite  strain.  Sun,  Ostien,  and  Salinger  (2013)  appar¬ 
ently  were  successful  in  extending  the  approach  of  White  and  Borja  (2008)  to  finite 
strain,  and  they  combined  it  with  an  assumed  enhanced  strain  approach  for  the  sol¬ 
id  skeleton  deformation  gradient  (for  near  incompressibility  of  the  solid  skeleton), 
which  is  an  additional  step  we  did  not  take  in  the  paper.  For  our  purposes,  the 
Brezzi  and  Pitkaranta  (1984)  approach  appears  effective,  and  we  will  use  it  in  the 
numerical  examples. 

After  applying  the  method  of  weighted  residuals  to  Eq.(14)  (Hughes,  1987),  substi¬ 
tuting  Darcy’s  law  from  Eq.(41),  dropping  the  (#)s  or  («)s  designation,  we  obtain 
the  variational  equation  of  the  mixture  balance  of  mass  in  the  reference  configura¬ 
tion  (of  the  solid  skeleton)  as, 


_  ovs?INT  i  os&INT  ,  n/pINT  ,  1  NT  ,  -^stab  ovs?EXT  r\ 

■7t  —  \  H-  e>^2  T"  ^3  \  e>^4  \  si  —  7/C^  —  U 

^NT=L"k 


\f))dV 


ANT 


=  [  %  p-i  fri  p-\ jjy 

Jm)dxln  dxK  Ki 


/}//7 

v  ty ?o  ”  — i  “  — 1\ 

<NT=l  ^T'b^-sdMV 

0  clXi 

JTstab=  /  ap-Fp^Fpj 
M  dx,  h  dxK  Kl 


1 JdV 


■7fXT  =  rjQfdA 


(47) 


where  a  is  the  stabilization  parameter,  and  <2f  is  the  normal  component  of  the  pore 
fluid  flux  (positive  inward)  across  the  boundary  Tq.  In  Truty  and  Zimmermann 
(2006),  an  analysis  was  conducted  to  determine  estimates  of  a  for  3D  problems, 
based  on  permeability,  solid  skeleton  stiffness,  pore  fluid  unit  weight,  time  step,  and 
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finite  element  size,  but  since  various  assumptions  were  made  on  other  factors,  we 
ended  up  estimating  a  by  a  trial-and-error  approach.  Further  investigation  may  be 
needed  to  more  efficiently  estimate  a  in  the  context  of  dynamic  loading  of  biphasic 
soft  tissues  at  finite  strain.  We  note  that  although  we  have  implemented  the 
term,  it  makes  little  difference  in  the  results,  and  then  may  be  left  out  of  future 
simulations. 

3.2  Element 

For  the  mixed  formulation,  we  use  a  Q27P8  hexahedral  element  as  shown  in  Fig.4 
(Regueiro  and  Ebrahimi,  2010).  The  stabilized  Q8P8  hexahedral  element  just  uses 
the  vertex  nodes  1-8. 


node  27  at  =  77  =  £  =  0 


Figure  4:  Q27P8:  27  nodes  for  tri-quadratic  interpolation  of  solid  skeleton  dis¬ 
placement  uh\  and  8  nodes  for  tri-linear  interpolation  of  pore  fluid  pressure  pf . 
Q8P8:  8  nodes  for  tri-linear  interpolation  of  solid  skeleton  displacement  uh\  and  8 
nodes  for  tri-linear  interpolation  of  pore  fluid  pressure  pf . 
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4  Numerical  examples 

Three  numerical  examples  are  presented  to  demonstrate  the  Q27P8  and  Q8P8  el¬ 
ements  (no  stabilization,  and  stabilized).  The  stable  Q27P8  formulation  is  imple¬ 
mented  in  the  FSSolidFluidMixT  class  in  Tahoe,  and  the  stabilized  Q8P8  ele¬ 
ment  is  implemented  in  the  FSSolidFluidMixQ8P8T  class  in  Tahoe  (stabilization 
is  turned  off  by  setting  a  —  0). 

4.1  Uniaxial  strain  compression  using  3D  FE 

We  start  with  a  uniaxial  strain  in  compression  example,  with  undrained  boundaries 
at  all  sides  except  the  top  boundary  (see  Fig. 5).  The  meshes  are  composed  of  ten 
Q27P8  and  ten  Q8P8  elements.  The  parameters  used  are  shown  in  Table  1. 


Table  1:  Parameters  for  lung  parenchyma  examples. 


source 

parameter 

units 

value 

Lande  and  Mitzner  (2006)  (air  through  parenchyma) 

ka 

m3.s/kg 

1  x  10“5 

Lande  and  Mitzner  (2006) 

«0 

0.99 

Lande  and  Mitzner  (2006)  (alveolar) 

^sK 

kg/m3 

1000 

Lande  and  Mitzner  (2006) 

V 

0.3 

Levental,  Georges,  and  Janmey  (2006)  (guinea  pig) 

E 

Pa 

5000 

this  paper 

Pa 

80,  1000,  10,000 

To  calculate  the  hydraulic  conductivity  kw  for  water,  we  consider  the  intrinsic  per¬ 
meability  calculated  from  the  hydraulic  conductivity  for  air  (Lande  and  Mitzner, 
2006 )x  =  kar\a  =  (1  x  10"5  m2/Pa.s)(1.83  x  10"5  Pa.s)  =  1.83  x  10"10  m2.  Then, 
we  use  the  viscosity  of  water  to  calculate  kw  =  >c/r\w=  1.83  x  10-10  m2/l  x  10-3 
Pa.s  =  1.83  x  10-7  m2/Pa.s.  The  bulk  moduli  of  air  and  water  are  Ka  —  1  x  105  Pa 
and  Kw  =  2.2  x  109  Pa,  respectively,  at  20°C. 

Referring  to  Fig. 5,  we  will  consider  three  types  of  loadings  and  solutions:  (1) 
drained  (quasi-static,  pore  fluid  pressure  pf  =  0,  tG  =  l,000Pa),  (2)  consolidat¬ 
ing  (coupled  pore  fluid  flow  and  solid  skeleton  deformation,  t°  =  l,000Pa),  and 
(3)  dynamic  impulse  loading  (coupled  pore  fluid  flow  and  solid  skeleton  deforma¬ 
tion  with  inertia  terms,  t°  =  10,000Pa).  For  each  loading  case,  we  keep  the  solid 
skeleton  parameters  the  same  for  the  lung  tissue,  and  then  assume  the  saturating 
fluid  is  either  (a)  air,  or  (b)  water,  with  parameters  given  in  Table  1 . 

For  pore  air,  the  solid  skeleton  vertical  displacement  upon  traction  loading  is  shown 
in  Fig.6(a).  It  can  be  seen  that  given  the  low  viscosity  of  the  pore  air,  the  drained 
and  consolidating  solutions  are  nearly  the  same,  meaning  the  pore  air  pressure  dis¬ 
sipates  as  the  column  is  loaded  over  0.1s,  as  illustrated  in  Fig. 6(b).  Nodes  8  and  91 
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drained 


consolidating  dynamic  impulse 


Figure  5:  Ten  hexahedral  element  column  mesh  for  the  three  analysis  cases:  (a) 
drained,  (b)  consolidating,  (c)  dynamic  impulse. 


are  the  locations  where  the  solutions  are  plotted  in  Fig.5.  In  Fig.6(b),  in  addition  to 
the  pore  fluid  pressure  pf  being  plotted  with  time,  the  mean  effective  Cauchy  stress 
p ',  and  effective  vertical  stress  o'zz  are  plotted.  Recall  that  ozz  and  p'  are  positive  in 
tension,  and  pf  is  positive  in  compression.  Gravity  is  ignored,  and  the  problem  is 
uniaxial  strain,  so  a'xx  and  c'yy  are  not  zero,  and  thus  p'  ^  Gzz. 

Next,  we  conduct  the  same  simulations,  but  switch  out  the  pore  fluid  from  air  to  wa¬ 
ter  properties  (viscosity,  compressibility).  The  solid  skeleton  vertical  displacement 
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pore  air  mid-length  corner  node  91  -  pore  air 


TIME  (sec)  TIME  (sec) 


(a)  (b) 

Figure  6:  (a)  Vertical  displacement,  and  (b)  Stress  for  drained  and  consolidating 
solutions  with  pore  air  flow. 


pore  water  mid-length  corner  node  91  -  pore  water 


(a)  (b) 

Figure  7:  (a)  Vertical  displacement,  and  (b)  Stress  for  drained  and  consolidating 
solutions  with  pore  water  flow. 


upon  traction  loading  is  shown  in  Fig. 7(a).  It  can  be  seen  that  given  the  higher  vis¬ 
cosity  of  the  pore  water,  the  drained  and  consolidating  solutions  are  now  different, 
as  it  takes  longer  for  the  pore  water  pressure  to  dissipate  with  time,  as  illustrated  in 
Fig. 7(b). 

Now,  to  investigate  the  significance  of  inertia  effects,  a  traction  load  of  t°  =  10, 000 
Pa  is  applied  over  0.1s  and  then  released  over  the  next  0.1s  for  the  pore  air  case 
(Fig.5),  and  loaded  over  0.01s  and  unloaded  over  0.01s  for  the  pore  water  case.  The 
effect  of  pore  fluid  (air  and  water)  is  investigated.  We  also  study  the  effect  of  global¬ 
ly  undrained  BC  as  compared  to  the  ends  of  the  column  being  drained.  In  Fig. 8(a), 


20  Copyright  ©  2014  Tech  Science  Press 


CMES  vol.98, ,  noA ,  pp.i-39,  2014 


we  plot  the  vertical  displacement  versus  time  for  nodes  8  and  91  for  drained  and 
undrained  BCs  for  the  coupled  dynamic  simulations  with  pore  air.  Recall  that  the 
bulk  modulus  of  air  at  ambient  temperature  is  relatively  low  (Ka  =  1  x  105Pa),  thus 
in  Fig.8(a)  there  is  compressive  vertical  displacement  even  for  the  undrained  BC 
case  (as  opposed  to  when  there  is  pore  water,  which  is  nearly  incompressible,  and 
thus  the  displacement  will  be  zero  for  the  undrained  BC  case  in  Fig.9(a)).  Also,  in 
Fig. 8(a),  note  that  for  the  case  with  drained  BC,  there  is  compressibility  of  pore  air 
and  also  relative  flow  of  pore  air,  reaching  a  displacement  near  0.06m,  which  for  an 
initial  column  length  of  0.1m,  is  nominally  =  0.06/0.1  =  60%  axial  strain,  which 
is  clearly  a  large  strain  (see  Fig.8(c)  for  actual  deformed  mesh  and  pf  contour).  In 
Fig. 8(b),  we  see  that  for  the  undrained  BC  case,  the  pore  fluid  pressure  pf  spikes 
up  to  nearly  10,000Pa  as  expected  (the  applied  traction  load),  while  because  of  the 
compressibility  of  air  and  relative  fluid  flow  even  during  globally  undrained  BCs, 
there  is  a  small  buildup  of  effective  stress. 


dynamic  impulse  -  pore  air 


P-f 

382.9 

300 

200 

:100 

6.008 


(a)  (b)  (c) 

Figure  8:  (a)  Vertical  displacement,  and  (b)  Stress  for  impulse  loading  solution  with 
pore  air  flow,  (c)  Actual  deformed  mesh  for  impulse  loading  with  pore  air  flow,  with 
drained  BCs.  Displacement  magnification  1  x .  Note  the  large  deformation.  Note 
the  original  aspect  ratio  in  Fig. 5  for  comparison. 


Next,  we  run  the  simulation  for  pore  water  with  10,000Pa  traction  over  0.01s,  with 
displacement  results  for  drained  and  undrained  BCs  shown  in  Fig. 9(a).  We  note 
now  what  we  expect  for  undrained  BC,  that  the  displacement  of  the  solid  skeleton 
is  zero  because  the  pore  water  is  nearly  incompressible  and  the  boundaries  are  im¬ 
permeable.  For  the  drained  BC,  we  observe  oscillation  of  the  solid  skeleton  vertical 
displacement  at  nodes  8  and  91,  which  eventually  damps  out  due  to  the  relative  flu¬ 
id  flow  effect.  Likewise,  in  Fig. 9(b),  we  observe  similar  expected  behavior,  where 
for  the  undrained  BC,  the  traction  is  completely  taken  up  by  the  pore  fluid  pressure 
Pf  ( tG  =  10,000Pa),  while  the  effective  stress  is  zero.  For  the  drained  BC,  pore 
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TIME  (sec)  TIME  (sec) 


(a)  (b) 

Figure  9:  (a)  Vertical  displacement,  and  (b)  Stress  for  impulse  loading  solution  with 
pore  water  flow. 


fluid  pressure  pf  builds  up  and  oscillates,  until  it  damps  out  to  near  zero  at  about 
0.2sec. 

Lastly,  we  compare  the  impulse  solution  as  if  the  material  were  completely  solid 
(no  relative  fluid  flow;  locally  drained),  with  results  shown  in  Fig.  10.  Since  the  real 
mass  densities  of  solid  and  fluid  (water)  are  the  same  for  the  lung  alveolar  tissue 
( psR  —  pfR  —  1000kg/m3),  the  initial  total  mass  density  is  also  po  =  1000kg/m3. 
We  can  compare  the  axial  wave  speed  estimated  from  the  curve  in  Fig.  10(a)  to  the 
small  strain  theoretical  solution  as  follows: 


(curve)  Vaxial  ft! 
(theory)  vaxiai  = 


2(0. lm) 

0.085sec 


=  2  Am / sec 


/A  +  2p  /(2885Ptf)  +  2(1923Ptf) 


1000 kg/m? 


—  2.6m/ sec 


(48) 

(49) 


Figure  10(b)  demonstrates  the  stress  solution  versus  time,  showing  that  for  the  sol¬ 
id  dynamic  solution  (no  pore  fluid  coupling),  we  get  the  expected  oscillatory  re¬ 
sponse,  with  little  algorithmic  damping  for  the  Newmark  time  integration  method 
used,  with  integration  parameters  /3  =  0.3025,  y  =  0.6,  unconditionally  stable,  with 
high  frequency  dissipation  (pg534  of  Hughes  (1987)).  Note  that  there  is  a  clear  dif¬ 
ference  between  solid  deformation  analysis  (locally  drained  material  points,  black 
curves)  and  coupled  pore  fluid  flow  with  solid  skeleton  deformation  analysis,  for 
both  undrained  BC  (red  curves)  and  drained  BC  (blue  curves).  The  drained  BC 
(blue  curve)  would  be  closest  to  the  actual  experimental  condition  in  the  lab,  where¬ 
as  the  undrained  BC  (red  curves)  could  be  replicated  if  the  lung  parenchyma  tis¬ 
sue  were  completely  sealed  along  its  boundaries  with  an  impermeable  membrane. 
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Clearly,  the  solid  analysis  (black  curves)  is  not  a  reasonable  assumption  for  the  tran¬ 
sient  behavior  of  a  soft  biological  tissue  treated  as  a  biphasic  (solid-fluid)  mixture 
at  finite  strain. 


TIME  (sec)  TIME  (sec) 


(a)  (b) 

Figure  10:  (a)  Vertical  displacement,  and  (b)  Stress  for  impulse  loading  solution 
with  pore  water  flow,  compared  to  solid  (locally  drained)  case. 


4.2  Uniaxial  strain  compression ,  with  stabilization 

We  now  re-run  the  simulations  for  pore  water  consolidation  and  impulse  loading, 
but  decrease  the  permeability  to  k  =  1  x  10“ 10  m2/Pa.s  to  test  the  stabilization  term 
in  Eq.(47)  for  a  =  1  x  10~6m3s2/kg.  Results  are  presented  in  Fig.  11  for  pore 
fluid  pressure  contours,  and  Fig.  12  for  pore  fluid  pressure  and  vertical  displacement 
versus  time.  We  see  in  the  contour  plots  the  stable  solution  for  Q27P8  element,  with 
near  uniform  pore  fluid  pressure  of  lOOOPa  in  the  middle  of  the  mesh  in  Fig.  11  (a), 
and  the  unstable  solution  with  the  Q8P8  element,  showing  oscillatory  pore  fluid 
pressure  along  the  depth  of  the  mesh  in  Fig.  11(b).  In  Fig.  11(c),  the  solution  is 
stabilized.  The  instability  is  evident  in  Fig.  12  too,  showing  apparent  consolidation, 
when  the  near  undrained  condition  caused  by  low  permeability  should  show  slow 
decrease  in  pore  fluid  pressure  pf  and  displacement  dz,  as  is  the  case  for  the  Q27P8 
element  and  the  stabilized  Q8P8  element. 

4.3  Pressure  impulse  loading  of  slab  of  lung  parenchyma 

The  second  numerical  example  considers  a  Neo-Hookean  isotropic  poroelastic  slab 
loaded  with  a  pressure  pulse  to  observe  transient  response  of  effective  stress  and 
pore  fluid  pressure  (Fig.  13).  The  slab  is  initially  assumed  to  be  saturated  with  water. 
The  top  face  of  the  slab  is  assumed  drained,  and  the  other  faces  are  impermeable. 
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(a)  (b)  (c) 


Figure  1 1 :  Contours  of  pore  fluid  pressure  pf  for  the  3  cases,  at  second  time  step  t  = 
0.1  sec:  (a)  Q27P8  element  showing  stable,  smooth  pf  contour;  (b)  Q8P8  element 
solution  showing  unstable,  oscillatory  pf  contour,  for  non- stabilized  solution;  (c) 
Q8P8  element  solution  with  stabilization  term. 


The  dynamic  loading  of  the  slab  as  biphasic  soft  tissue  is  shown  in  Fig.  13.  The 
load-time  schedule  functions  for  drained,  consolidating  (same  as  load  and  hold  for 
dynamic  simulation),  and  dynamic  load  and  release  (impulse)  are  also  illustrated  in 
Fig.  13.  For  constitutive  parameters,  we  use  the  same  as  in  Table  1.  For  the  Q27P8 
element,  50  and  200  element  meshes  are  considered.  For  the  Q8P8  element,  a  1350 
element  mesh  is  considered. 

First,  consider  the  drained  and  consolidating  deformed  meshes  in  Fig.  14.  This 
gives  a  general  deformation  pattern  for  the  deformed  meshes.  Note  nearly  the  same 
deformation  is  obtained  with  the  3  meshes.  Next,  consider  the  deformed  meshes 
for  dynamic  simulation  for  hold  and  impulse  loadings  in  Fig.  15  for  pore  water,  for 
the  50  Q27P8  element  mesh.  The  hold  loading  generates  more  deformation  and 
also  higher  peak  pore  fluid  pressure  than  the  impulse  loading,  as  shown  in  Fig.  16. 
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Figure  12:  (a)  Pore  water  pressure  pf  versus  time  at  middle  section  of  the  mesh, 
comparing  result  of  Q27P8  mesh  (stable,  blue  lines)  and  Q8P8  meshes  (non- 
stabilized  (red)  and  stabilized  (green)),  (b)  Displacement  in  z  direction  dz  versus 
time  at  top  and  middle  sections  of  the  mesh,  comparing  result  of  Q27P8  mesh  (sta¬ 
ble)  and  Q8P8  meshes  (non- stabilized  and  stabilized). 


In  Fig.  16,  we  compare  the  pore  fluid  pressure  pf  at  the  corner  node  for  the  various 
analysis  types.  We  see  that  for  the  dynamic  hold  loading,  pf  is  much  larger  than  for 
dynamic  impulse  loading.  For  all  drained  solutions,  pf  =  0,  and  for  the  consolidat¬ 
ing  solution,  there  is  a  small  negative  pf  before  immediate  consolidation  (the  solid 
skeleton  at  the  bottom  impermeable  face  is  in  tension).  The  200  Q27P8  element 
mesh  gives  almost  the  exact  same  result  as  in  Fig.  16(a)  (and  is  thus  not  shown), 
and  the  non- stabilized  1350  Q8P8  element  mesh  gives  nearly  the  same  result  as  in 
Fig.  16(a),  with  result  shown  in  Fig.  16(b).  Thus,  for  some  soft  biological  tissues 
with  higher  permeability  (such  as  the  lung  parenchyma,  not  the  case  for  the  low 
permeability  vertebral  disc),  stable  results  can  be  achieved  with  the  non- stabilized 
Q8P8  element  implementation.  In  Fig.  17,  we  compare  displacement  dy  at  the  cor¬ 
ner  node  for  the  various  analyses,  plotting  over  0.1  sec  and  lsec.  We  see  that  the 
excited  displacement  frequencies  are  different  for  the  dynamic  hold  and  impulse 
simulations,  the  hold  one  having  a  higher  frequency  partly  due  to  stiffening  by  en¬ 
abling  geometric  nonlinearity.  We  see  in  Fig.  17(c),  the  1350  Q8P8  element  mesh 
result  is  nearly  the  same  as  the  50  Q27P8  element  mesh  result  in  Fig.  17(a).  Given 
the  small  pore  pressure  pf  generated  during  impulse  loading  (red  curves  in  Fig.  16), 
there  is  no  noticeable  difference  between  “dynamic  impulse”  and  “dynamic  im¬ 
pulse  drained”  displacement  time  histories  in  Fig.  17,  whereas  for  “dynamic  hold” 
and  “dynamic  hold  drained,”  there  is  a  small  difference  in  displacement  amplitude. 
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drained  consolidating  dynamic  impulse 

to  =  1ms,  t\  =  2ms 


Figure  13:  (top)  Schematic  of  pressure  loading  of  cantilevered  poroelastic  slab, 
(bottom)  Various  load  schedule  functions  applied. 
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Figure  14:  The  deformed  mesh  for  drained  loading,  and  end  of  consolidation  (no 
inertia  terms)  for  pore  water.  Displacement  units  in  meters,  with  displacement 
magnification  lx.  (a)  50  Q27P8  elements,  (b)  200  Q27P8  elements,  (c)  1350  Q8P8 
elements. 
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(a)  (b) 

Figure  15:  The  deformed  meshes  for  dynamic  simulation  with  pore  water  for  hold 
(a)  and  impulse  (b)  loadings  (with  inertia  terms),  pf  units  in  Pa. 
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Figure  16:  Plot  of  pore  water  pressure  pf  versus  time  for  the  various  analyses, 
(a)  50  Q27P8  element  mesh  results,  (b)  1350  Q8P8  element  mesh  results  without 
stabilization. 
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Figure  17:  Plot  of  displacement  dy  versus  time  for  the  various  analyses  ((a)  up  to 
O.lsec,  (b)  up  to  lsec  for  50  Q27P8  element  mesh)  for  pore  water,  (c)  1350  Q8P8 
element  mesh  results  without  stabilization. 


Next,  we  consider  the  deformed  meshes  for  dynamic  simulation  for  impulse  and 
hold  loadings  in  Fig.  18  for  pore  air.  In  this  case  for  pore  air,  the  peak  pore  air  pres¬ 
sures  are  not  noticeably  different  between  dynamic  hold  and  impulse  loadings,  as 
also  shown  in  Fig.  19(a)  for  pf  versus  time.  In  Fig.  19(b),  we  compare  displacemen- 
t  dy  at  the  corner  node,  plotting  over  O.lsec  for  pore  air.  We  see  that  the  dynamic 
hold  and  impulse  simulations  have  nearly  the  same  frequency,  while  there  is  a  small 
poromechanical  effect  noticeable  for  the  dynamic  impulse  simulation  (difference  in 
curves  after  0.05sec).  For  the  Q8P8  element  mesh,  the  results  are  similar  as  shown 
in  Fig.20. 
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(a)  (b) 

Figure  18:  The  deformed  meshes  for  dynamic  simulation  with  pore  air  for  hold  (a) 
and  impulse  (b)  loadings  (with  inertia  terms),  pf  units  in  Pa. 
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Figure  19:  50  Q27P8  element  mesh  results:  (a)  Plot  of  pore  air  pressure  pf  versus 
time  for  the  various  analyses,  (b)  Plot  of  displacement  dy  versus  time  for  the  various 
analyses. 
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Figure  20:  Non- stabilized  1350  Q8P8  element  mesh  results:  (a)  Plot  of  pore  air 
pressure  pf  versus  time  for  the  various  analyses,  (b)  Plot  of  displacement  dy  versus 
time  for  the  various  analyses. 
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In  summary,  for  longer  duration  loadings  (in  this  case,  up  to  0. 1  sec),  the  pore  water 
pressure  can  build  up  and  thus  affect  the  solid  skeleton  deformation,  and  in  turn  the 
effective  stresses.  When  the  pore  water  is  substituted  with  pore  air,  differences  in 
results  are  observed  between  the  two  pore  fluids.  For  pore  water  (higher  mixture 
mass  density),  different  loadings  (dynamic  impulse,  dynamic  hold)  trigger  different 
modes  of  vibration  at  different  frequencies.  On  the  other  hand,  for  pore  air  (lower 
mixture  mass  density),  different  loadings  (dynamic  impulse,  dynamic  hold)  trigger 
the  same  mode  of  vibration  (with  slightly  different  frequencies  due  to  geometric 
nonlinearities).  For  pore  air,  the  overall  mixture  is  less  dense,  thus  the  frequency 
of  free  or  forced  vibration  is  higher  than  when  saturated  with  pore  water;  the  per¬ 
meability  is  55  x  higher  than  for  pore  water,  thus  there  is  little  build  up  of  pore  air 
pressure  during  dynamic  loading  of  the  thin  slab. 

Also,  it  was  noted  that  there  is  no  distinguishable  difference  between  a  50  element 
and  200  element  mesh  for  Q27P8  implementation,  and  likewise  for  a  1350  element 
Q8P8  mesh  without  stabilization.  Because  of  the  higher  permeability  of  the  lung 
parenchyma,  the  stabilizing  term  ^ stab  is  not  needed.  We  will  see  in  the  next 
example  for  the  vertebral  disc  that  has  a  much  smaller  permeability,  the  stabilizing 
term  is  needed  to  obtain  meaningful,  stable  results. 

4.4  Vertebral  disc  compression  loading 

Vertebral  disc  compression  using  the  Q27P8  element,  and  Q8P8  element  with  and 
without  stabilization,  is  simulated  with  parameters  in  Table  2  (taken  from  Tables  1 
and  2  of  Williams,  Natarajan,  and  Andersson  (2007))  and  coarsest  mesh  in  Fig.21; 
two  finer  Q8P8  meshes  are  shown  in  Fig. 22.  The  geometry  of  the  lumbar  spine 
disc  was  simplified  for  the  simulations.  The  annulus  fibrosus  was  modeled  as  an 
ellipse  with  a  major  axis  of  0.02646  m  and  a  minor  axis  of  0.02016m  with  a  height 
of  0.00945  m  (minus  the  nucleus  volume).  The  nucleus  pulposus  was  modeled  as 
an  ellipse  with  a  major  axis  of  0.02016  m  and  a  minor  axis  of  0.01386  m  with  a 
height  of  0.00945  m.  The  permeabilities  k  are  averaged  from  the  values  reported 
in  different  directions  in  Williams,  Natarajan,  and  Andersson  (2007),  assuming 
anisotropy  in  that  paper;  we  assume  isotropic  permeability  in  this  paper,  which  can 
be  generalized  for  anisotropy  in  the  future.  A  downward  traction  load  of  0.5  MPa 
(Ferguson,  Ito,  and  Nolte,  2004)  is  applied,  with  similar  time  histories  as  shown 
in  Fig.5  for  drained,  consolidating,  and  dynamic  impulse.  We  do  not  verify  that 
the  material  properties  in  Williams,  Natarajan,  and  Andersson  (2007)  for  annulus 
and  nucleus  of  the  disc,  nor  the  compressive  loading  of  0.5MPa  by  Ferguson,  Ito, 
and  Nolte  (2004)  are  physiologically  correct  or  not.  The  purpose  of  the  paper  is 
to  investigate  the  difference  in  results  in  assuming  a  drained  versus  consolidating 
versus  dynamic  (with  inertia  terms)  coupled  poromechanical  analysis.  All  analyses 
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use  stabilization  parameter  a  =  1  x  10  6m3s2 /kg. 


Table  2:  Parameters  for  vertebral  disc  compression  (Williams,  Natarajan,  and  An- 
dersson,  2007). 


tissue 

parameter 

units 

value 

nucleus 

k 

m3.s/kg 

1.8x  10-15 

nucleus 

4 

0.83 

nucleus 

pSR 

kg/m3 

1000 

nucleus 

A 

Pa 

3.1  xlO6 

nucleus 

M 

Pa 

0.345  xlO6 

annulus 

k 

m3.s/kg 

1.7xl0“15 

annulus 

n0 

0.78 

annulus 

pSR 

kg/m3 

1060 

annulus 

A 

Pa 

3.6xl06 

annulus 

Pa 

0.893  xlO6 

Figure  21:  Coarsest  mesh  (432  elements)  used  in  analyzing  vertebral  disc  compres¬ 
sion,  showing  applied  effective  traction  t°'  and  drained  pf  =  0  boundary  at  top,  and 
Cases  1  and  2  BCs  on  the  curved  surface  of  the  annulus.  The  nucleus  (blue)  and 
annulus  (red)  sections  are  modeled  as  distinct  materials.  The  symmetry  planes  are 
impermeable  with  fixed  normal  displacements. 
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4.4.1  Vertebral  disc  compression  loading  with  lateral  expansion  and  drainage  - 
Case  1 

For  the  Case  1  BCs  ( t° '  =  0,  pf  =  0  on  curved  annulus  surface)  shown  in  Fig.21, 
we  present  the  time  histories  of  displacement  dz  (vertical,  in  direction  of  traction 
loading)  and  dx  (lateral,  along  long  axis),  and  pore  fluid  pressure  /?f,  in  Fig. 25  for 
the  three  nodes  shown  in  Fig.23.  The  time  histories  in  Fig.25  provide  a  comparison 
between  drained,  consolidating,  and  dynamic  impulse  results,  as  well  as  various 
meshes  (Q27P8  and  Q8P8  element  types).  Unless  otherwise  indicated,  it  is  as¬ 
sumed  the  result  is  that  of  the  coarsest  mesh  (432  elements)  with  Q27P8  element 
type.  In  Fig.24,  we  see  contour  plots  of  pf  with  deformed  meshes  for  the  consol¬ 
idating  analysis.  Note  from  Fig.25,  that  upon  end  of  traction  loading  at  0.1s,  the 
displacements  dz  and  dx  and  pore  fluid  pressure  pf  remain  nearly  constant,  where 
in  fact  there  is  a  slow  consolidation  process  ongoing  (note  the  small  permeabil¬ 
ity  values  in  Table  2  for  nucleus  and  annulus).  Another  observation  is  that  the 
Q27P8  mesh  provides  a  higher  pf ,  yet  with  nearly  the  same  displacement  as  the 
stabilized  Q8P8  meshes  (see  Fig.25).  It  is  likely  that  in  the  stabilized  Q8P8  mesh 
results,  there  is  added  diffusion  of  pore  fluid  leading  to  lower  pf  values  as  observed 
in  Fig.24.  Mesh  refinement  on  the  Q27P8  mesh  is  required  for  further  study,  but 
given  the  expense  of  this  element  type,  and  the  lack  of  parallel  execution  for  the 
Q27P8  currently,  we  do  not  present  finer  mesh  results  for  the  Q27P8  element.  The 
Q8P8  element  can  be  executed  in  parallel  computation.  For  the  dynamic  impulse 
analysis,  given  the  three-dimensional  nature  of  loading  and  pore  fluid  flow,  the  disc 
essentially  damps  out  the  impulse  wave  completely,  with  small  oscillation  after  re¬ 
lease  of  the  impulse  at  0.02s,  as  shown  in  Fig.25(a),(b)  displacements.  The  drained 


Figure  22:  Two  finer  Q8P8  meshes  used  in  the  disc  compression  analysis,  showing 
nodes  where  analysis  results  are  compared.  Left  mesh  has  3,465  elements,  and  the 
right  mesh  has  16,875  elements. 
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results  (no  consolidation,  and  no  inertia  terms)  are  provided  as  a  benchmark  for 
deformation  comparison. 


d_  Magnitude 


Figure  23:  Q27P8  deformed  mesh  for  drained  analysis  of  Case  1  BC  ( t° '  =  0,  pf  = 
0).  Observe  the  undeformed  mesh  for  comparison.  Displacement  magnification 
lx. 


pore_press  pore_press 


(b)  (c) 


Figure  24:  Deformed  meshes  (displacement  magnification  lx),  with  contour  of 


pore  fluid  pressure  pf ,  for  the  following  meshes:  (a)  coarse  Q27P8  mesh  (432  el¬ 
ements),  (b)  stabilized  coarse  Q8P8  mesh  (432  elements),  (c)  stabilized  fine  Q8P8 


mesh  (3465  elements). 
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drained-4307 
d  rai  ned-hex8coarse-662 
drained-hex8fine-441 2 
drained-hex8finer-1 8962 
consol-4307 
consol-hex8coarse-662 
consol-hex8fine-441 2 
impulse-4307 
i  mpulse-hex8coarse-662 
impulse-hex8fine-441 2 


0.1 

time  (sec) 


(a) 


(b) 


time  (sec) 

(c) 


time  (sec) 

(d) 


Figure  25:  (a)  Vertical  displacement  dz.  (b)  Lateral  displacement  dx  along  long 
axis,  (c)  Pore  fluid  pressure  pf.  (d)  Pore  fluid  pressure  pf  without  Q27P8  result. 


d_  Magnitude 
0.0008973 
10.0008 
;0.0006 
0.0004 
0.0002 

0 


Figure  26:  Q27P8  deformed  mesh  for  drained  analysis  of  Case  2  BC  ( un  =  0,  Sw  = 
0).  Observe  the  undeformed  mesh  for  comparison.  Displacement  magnification 
lx. 
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(C)  (d) 

Figure  27:  Deformed  meshes  (displacement  magnification  lx),  with  contour  of 
pore  fluid  pressure  pf,  for  the  following  meshes:  (a)  coarse  Q27P8  mesh  (432  ele¬ 
ments),  (b)  non- stabilized  coarse  Q8P8  mesh  (432  elements),  (c)  stabilized  coarse 
Q8P8  mesh  (432  elements),  (d)  stabilized  fine  Q8P8  mesh  (3465  elements). 
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(a) 
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time  (sec) 


(b)  (c) 

Figure  28:  (a)  Vertical  displacement  dz.  (b)  Pore  fluid  pressure  pf.  (c)  Pore  fluid 
pressure  pf  without  Q27P8  result. 


4.4.2  Vertebral  disc  compression  loading  in  uniaxial  strain  -  Case  2 

For  the  Case  2  BCs  ( un  =  0,  Sw  =  0  on  curved  annulus  surface)  shown  in  Fig.21, 
we  present  the  time  histories  of  displacement  dz  and  dx ,  and  pore  fluid  pressure 
Pi ,  in  Fig. 28  for  the  same  three  nodes  shown  in  Fig.23,  with  drained  deformed 
mesh  solution  shown  in  Fig.26.  Here,  the  higher  pore  fluid  pressure  pf  generated 
in  the  Q27P8  element  mesh  indicates  a  near  locking  behavior  because  of  the  low 
permeability,  that  stiffens  the  overall  compressibility  of  the  disc  such  that  dz  is  less 
for  the  Q27P8  mesh  versus  the  Q8P8  meshes  with  stabilization  (see  Fig. 28(a)).  The 
two  stabilized  Q8P8  meshes  compare  well  with  each  other  as  seen  in  Fig. 28(a), (c). 
Also  note  for  non- stabilized  Q8P8  coarse  mesh,  the  classical  oscillating  pore  fluid 
pressure  is  observed  in  Fig. 27(b),  which  the  stabilized  term  alleviates  in  Fig. 27(c). 

Similar  to  Case  1,  there  is  only  a  small  oscillation  in  solid  skeleton  displacement 
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upon  release  of  the  impulse  traction,  implying  that  the  forced  solid  skeleton  com¬ 
pression  and  driven  pore  fluid  flow  almost  completely  damps  out  the  impulse  load. 

5  Conclusions 

The  paper  presented  finite  strain  three-dimensional  finite  element  analysis  (FEA) 
of  coupled  pore  fluid  flow  and  solid  skeleton  deformation  of  soft  biological  tissues 
at  finite  strain,  including  wave  propagation  (i.e.,  inertia  terms).  Compressibility  of 
the  pore  fluid  is  considered,  as  well  as  porosity  dependent  permeability.  The  FEA 
results  of  a  full  mixed  formulation  Q27P8  (triquadratic  in  displacement,  trilinear  in 
pore  fluid  pressure)  hexahedral  element  is  compared  to  an  equal  order  interpolation 
mixed  formulation  Q8P8  (trilinear  in  displacement,  trilinear  in  pore  fluid  pressure) 
hexahedral  element  with  stabilization.  It  was  shown  that  the  stabilization  method  of 
Brezzi  and  Pitkaranta  (1984)  is  effective  for  finite  strain  coupled  FEA  of  biphasic 
mixtures  of  soft  tissues  with  small  permeability  (such  as  the  vertebral  disc),  and 
is  not  needed  for  higher  permeability  tissues  (such  as  the  lung  parenchyma).  It  is 
useful  to  have  the  Q27P8  element  implementation  (even  though  computationally 
expensive),  as  a  benchmark  for  comparison  to  the  Q8P8  element  implementation 
(non- stabilized,  and  stabilized). 

Two  soft  biological  tissues  are  considered:  (1)  lung  parenchyma,  and  (2)  vertebral 
disc.  For  (1),  it  was  observed  that  switching  out  the  pore  fluid  between  water  and 
air  leads  to  different  dynamical  results,  where  the  less  viscous  and  lighter  (but  more 
compressible)  pore  air  can  nearly  damp  out  any  applied  impulse  wave,  but  also  lead 
to  different  excited  mode  shapes  as  a  result  of  the  lower  mass  density.  For  the  pore 
water  case,  some  waves  can  be  observed  after  application  of  an  impulse  traction 
load,  but  they  are  damped  quickly  due  to  the  dissipation  associated  with  relative 
pore  fluid  flow  and  solid  skeleton  deformation.  For  (2),  given  the  low  permeability 
of  the  disc,  and  assuming  saturation  by  pore  water,  application  of  traction  dynamic 
impulse  load  leading  to  large  compressive  solid  skeleton  deformation  (and  corre¬ 
sponding  relative  pore  fluid  flow)  does  not  allow  any  significant  waves  to  propagate 
in  time,  providing  a  near  fully  damped  response,  but  with  observable  change  in  pore 
fluid  pressure. 

It  is  clear  that  there  is  an  interplay  of  solid  skeleton  deformation  (and,  in  turn,  solid 
skeleton  effective  stress)  and  build-up  of  pore  fluid  pressure  over  time,  that  would 
not  be  properly  accounted  for  by  a  more  simple  finite  strain  viscoelastic  constitutive 
model,  without  treating  the  soft  tissue  as  a  biphasic  mixture.  The  effective  stress 
governs  the  constitutive  response  of  the  solid  skeleton,  which  is  influenced  by  the 
changing  pore  fluid  pressure  with  time.  Thus,  if  appropriate  anisotropic,  damage 
constitutive  models  are  to  be  considered  for  soft  biological  tissues  loaded  dynami¬ 
cally  (such  as  during  impact  or  blast  loading),  the  biphasic  nature  of  the  soft  tissue 


Soft  Biological  Tissues  Treated  as  Biphasic  Porous  Media 


37 


should  be  taken  into  account.  With  that  said,  the  paper  provides  a  glimpse  into  how 
such  biphasic  dynamic  (with  inertia  terms)  analyses  may  be  conducted,  benchmark¬ 
ing  against  drained  and  consolidating  FEA  results,  with  further  research  needed:  (i) 
consider  physiological  BCs  and  tissue  geometries  for  in-vivo  conditions;  (ii)  more 
study  of  the  stabilization  parameter  a  and  how  to  select  it  for  various  tissue  elastic 
moduli,  geometries,  etc.;  (iii)  anisotropic  permeability  and  solid  skeleton  consti¬ 
tutive  models  appropriate  for  fibrous  soft  biological  tissues;  (iv)  inclusion  of  pore 
fluid  acceleration  vector  different  than  solid  skeleton  acceleration  (i.e.,  ctf  as)  and 
modification  of  the  coupled  governing  equation  solution  algorithm  (see  Remark  1 
in  Sect.2.4);  and  (v)  extension  to  dynamic  explicit  analysis  capability  for  high  rate 
dynamic  loading. 
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